// eee.Sheffield.PZ.Math
//
// Copyright ?Ping Zou, 2007
// sg71.cherub@gmail.com

using System;
using System.Collections.Generic;
using System.Text;

namespace eee.Sheffield.PZ.Math
{
    struct table 
    {
        public int _n;
        public double _f;
        public long _i;

        public table(int n, double f, long i)
        {
            _n = n;
            _f = f;
            _i = i;
        }
    }
    
    #region GSL comments
    /* specfunc/gamma.c
     * 
     * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
     * 
     * This program is free software; you can redistribute it and/or modify
     * it under the terms of the GNU General Public License as published by
     * the Free Software Foundation; either version 2 of the License, or (at
     * your option) any later version.
     * 
     * This program is distributed in the hope that it will be useful, but
     * WITHOUT ANY WARRANTY; without even the implied warranty of
     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     * General Public License for more details.
     * 
     * You should have received a copy of the GNU General Public License
     * along with this program; if not, write to the Free Software
     * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
     */

    /* Author:  G. Jungman */
    #endregion
    public class Gamma
    {
        #region static parameters
        public static double LogRootTwoPi_ = 0.9189385332046727418;
        static int GSL_SF_FACT_NMAX = 170;
        static int GSL_SF_DOUBLEFACT_NMAX = 297;
        static int GSL_SF_GAMMA_XMAX = 171;
        static table[] fact_table = {
            new table ( 0,  1.0,     1L     ),
            new table ( 1,  1.0,     1L     ),
            new table ( 2,  2.0,     2L     ),
            new table ( 3,  6.0,     6L     ),
            new table ( 4,  24.0,    24L    ),
            new table ( 5,  120.0,   120L   ),
            new table ( 6,  720.0,   720L   ),
            new table ( 7,  5040.0,  5040L  ),
            new table ( 8,  40320.0, 40320L ),

            new table ( 9,  362880.0,     362880L    ),
            new table ( 10, 3628800.0,    3628800L   ),
            new table ( 11, 39916800.0,   39916800L  ),
            new table ( 12, 479001600.0,  479001600L ),

            new table ( 13, 6227020800.0,                               0 ),
            new table ( 14, 87178291200.0,                              0 ),
            new table ( 15, 1307674368000.0,                            0 ),
            new table ( 16, 20922789888000.0,                           0 ),
            new table ( 17, 355687428096000.0,                          0 ),
            new table ( 18, 6402373705728000.0,                         0 ),
            new table ( 19, 121645100408832000.0,                       0 ),
            new table ( 20, 2432902008176640000.0,                      0 ),
            new table ( 21, 51090942171709440000.0,                     0 ),
            new table ( 22, 1124000727777607680000.0,                   0 ),
            new table ( 23, 25852016738884976640000.0,                  0 ),
            new table ( 24, 620448401733239439360000.0,                 0 ),
            new table ( 25, 15511210043330985984000000.0,               0 ),
            new table ( 26, 403291461126605635584000000.0,              0 ),
            new table ( 27, 10888869450418352160768000000.0,            0 ),
            new table ( 28, 304888344611713860501504000000.0,           0 ),
            new table ( 29, 8841761993739701954543616000000.0,          0 ),
            new table ( 30, 265252859812191058636308480000000.0,        0 ),
            new table ( 31, 8222838654177922817725562880000000.0,       0 ),
            new table ( 32, 263130836933693530167218012160000000.0,     0 ),
            new table ( 33, 8683317618811886495518194401280000000.0,    0 ),
            new table ( 34, 2.95232799039604140847618609644e38,  0 ),
            new table ( 35, 1.03331479663861449296666513375e40,  0 ),
            new table ( 36, 3.71993326789901217467999448151e41,  0 ),
            new table ( 37, 1.37637530912263450463159795816e43,  0 ),
            new table ( 38, 5.23022617466601111760007224100e44,  0 ),
            new table ( 39, 2.03978820811974433586402817399e46,  0 ),
            new table ( 40, 8.15915283247897734345611269600e47,  0 ),
            new table ( 41, 3.34525266131638071081700620534e49,  0 ),
            new table ( 42, 1.40500611775287989854314260624e51,  0 ),
            new table ( 43, 6.04152630633738356373551320685e52,  0 ),
            new table ( 44, 2.65827157478844876804362581101e54,  0 ),
            new table ( 45, 1.19622220865480194561963161496e56,  0 ),
            new table ( 46, 5.50262215981208894985030542880e57,  0 ),
            new table ( 47, 2.58623241511168180642964355154e59,  0 ),
            new table ( 48, 1.24139155925360726708622890474e61,  0 ),
            new table ( 49, 6.08281864034267560872252163321e62,  0 ),
            new table ( 50, 3.04140932017133780436126081661e64,  0 ),
            new table ( 51, 1.55111875328738228022424301647e66,  0 ),
            new table ( 52, 8.06581751709438785716606368564e67,  0 ),
            new table ( 53, 4.27488328406002556429801375339e69,  0 ),
            new table ( 54, 2.30843697339241380472092742683e71,  0 ),
            new table ( 55, 1.26964033536582759259651008476e73,  0 ),
            new table ( 56, 7.10998587804863451854045647464e74,  0 ),
            new table ( 57, 4.05269195048772167556806019054e76,  0 ),
            new table ( 58, 2.35056133128287857182947491052e78,  0 ),
            new table ( 59, 1.38683118545689835737939019720e80,  0 ),
            new table ( 60, 8.32098711274139014427634118320e81,  0 ),
            new table ( 61, 5.07580213877224798800856812177e83,  0 ),
            new table ( 62, 3.14699732603879375256531223550e85,  0 ),
            new table ( 63, 1.982608315404440064116146708360e87,  0 ),
            new table ( 64, 1.268869321858841641034333893350e89,  0 ),
            new table ( 65, 8.247650592082470666723170306800e90,  0 ),
            new table ( 66, 5.443449390774430640037292402480e92,  0 ),
            new table ( 67, 3.647111091818868528824985909660e94,  0 ),
            new table ( 68, 2.480035542436830599600990418570e96,  0 ),
            new table ( 69, 1.711224524281413113724683388810e98,  0 ),
            new table ( 70, 1.197857166996989179607278372170e100,  0 ),
            new table ( 71, 8.504785885678623175211676442400e101,  0 ),
            new table ( 72, 6.123445837688608686152407038530e103,  0 ),
            new table ( 73, 4.470115461512684340891257138130e105,  0 ),
            new table ( 74, 3.307885441519386412259530282210e107,  0 ),
            new table ( 75, 2.480914081139539809194647711660e109,  0 ),
            new table ( 76, 1.885494701666050254987932260860e111,  0 ),
            new table ( 77, 1.451830920282858696340707840860e113,  0 ),
            new table ( 78, 1.132428117820629783145752115870e115,  0 ),
            new table ( 79, 8.946182130782975286851441715400e116,  0 ),
            new table ( 80, 7.156945704626380229481153372320e118,  0 ),
            new table ( 81, 5.797126020747367985879734231580e120,  0 ),
            new table ( 82, 4.753643337012841748421382069890e122,  0 ),
            new table ( 83, 3.945523969720658651189747118010e124,  0 ),
            new table ( 84, 3.314240134565353266999387579130e126,  0 ),
            new table ( 85, 2.817104114380550276949479442260e128,  0 ),
            new table ( 86, 2.422709538367273238176552320340e130,  0 ),
            new table ( 87, 2.107757298379527717213600518700e132,  0 ),
            new table ( 88, 1.854826422573984391147968456460e134,  0 ),
            new table ( 89, 1.650795516090846108121691926250e136,  0 ),
            new table ( 90, 1.485715964481761497309522733620e138,  0 ),
            new table ( 91, 1.352001527678402962551665687590e140,  0 ),
            new table ( 92, 1.243841405464130725547532432590e142,  0 ),
            new table ( 93, 1.156772507081641574759205162310e144,  0 ),
            new table ( 94, 1.087366156656743080273652852570e146,  0 ),
            new table ( 95, 1.032997848823905926259970209940e148,  0 ),
            new table ( 96, 9.916779348709496892095714015400e149,  0 ),
            new table ( 97, 9.619275968248211985332842594960e151,  0 ),
            new table ( 98, 9.426890448883247745626185743100e153,  0 ),
            new table ( 99, 9.332621544394415268169923885600e155,  0 ),
            new table ( 100, 9.33262154439441526816992388563e157,  0 ),
            new table ( 101, 9.42594775983835942085162312450e159,  0 ),
            new table ( 102, 9.61446671503512660926865558700e161,  0 ),
            new table ( 103, 9.90290071648618040754671525458e163,  0 ),
            new table ( 104, 1.02990167451456276238485838648e166,  0 ),
            new table ( 105, 1.08139675824029090050410130580e168,  0 ),
            new table ( 106, 1.146280563734708354534347384148e170,  0 ),
            new table ( 107, 1.226520203196137939351751701040e172,  0 ),
            new table ( 108, 1.324641819451828974499891837120e174,  0 ),
            new table ( 109, 1.443859583202493582204882102460e176,  0 ),
            new table ( 110, 1.588245541522742940425370312710e178,  0 ),
            new table ( 111, 1.762952551090244663872161047110e180,  0 ),
            new table ( 112, 1.974506857221074023536820372760e182,  0 ),
            new table ( 113, 2.231192748659813646596607021220e184,  0 ),
            new table ( 114, 2.543559733472187557120132004190e186,  0 ),
            new table ( 115, 2.925093693493015690688151804820e188,  0 ),
            new table ( 116, 3.393108684451898201198256093590e190,  0 ),
            new table ( 117, 3.96993716080872089540195962950e192,  0 ),
            new table ( 118, 4.68452584975429065657431236281e194,  0 ),
            new table ( 119, 5.57458576120760588132343171174e196,  0 ),
            new table ( 120, 6.68950291344912705758811805409e198,  0 ),
            new table ( 121, 8.09429852527344373968162284545e200,  0 ),
            new table ( 122, 9.87504420083360136241157987140e202,  0 ),
            new table ( 123, 1.21463043670253296757662432419e205,  0 ),
            new table ( 124, 1.50614174151114087979501416199e207,  0 ),
            new table ( 125, 1.88267717688892609974376770249e209,  0 ),
            new table ( 126, 2.37217324288004688567714730514e211,  0 ),
            new table ( 127, 3.01266001845765954480997707753e213,  0 ),
            new table ( 128, 3.85620482362580421735677065923e215,  0 ),
            new table ( 129, 4.97450422247728744039023415041e217,  0 ),
            new table ( 130, 6.46685548922047367250730439554e219,  0 ),
            new table ( 131, 8.47158069087882051098456875820e221,  0 ),
            new table ( 132, 1.11824865119600430744996307608e224,  0 ),
            new table ( 133, 1.48727070609068572890845089118e226,  0 ),
            new table ( 134, 1.99294274616151887673732419418e228,  0 ),
            new table ( 135, 2.69047270731805048359538766215e230,  0 ),
            new table ( 136, 3.65904288195254865768972722052e232,  0 ),
            new table ( 137, 5.01288874827499166103492629211e234,  0 ),
            new table ( 138, 6.91778647261948849222819828311e236,  0 ),
            new table ( 139, 9.61572319694108900419719561353e238,  0 ),
            new table ( 140, 1.34620124757175246058760738589e241,  0 ),
            new table ( 141, 1.89814375907617096942852641411e243,  0 ),
            new table ( 142, 2.69536413788816277658850750804e245,  0 ),
            new table ( 143, 3.85437071718007277052156573649e247,  0 ),
            new table ( 144, 5.55029383273930478955105466055e249,  0 ),
            new table ( 145, 8.04792605747199194484902925780e251,  0 ),
            new table ( 146, 1.17499720439091082394795827164e254,  0 ),
            new table ( 147, 1.72724589045463891120349865931e256,  0 ),
            new table ( 148, 2.55632391787286558858117801578e258,  0 ),
            new table ( 149, 3.80892263763056972698595524351e260,  0 ),
            new table ( 150, 5.71338395644585459047893286526e262,  0 ),
            new table ( 151, 8.62720977423324043162318862650e264,  0 ),
            new table ( 152, 1.31133588568345254560672467123e267,  0 ),
            new table ( 153, 2.00634390509568239477828874699e269,  0 ),
            new table ( 154, 3.08976961384735088795856467036e271,  0 ),
            new table ( 155, 4.78914290146339387633577523906e273,  0 ),
            new table ( 156, 7.47106292628289444708380937294e275,  0 ),
            new table ( 157, 1.17295687942641442819215807155e278,  0 ),
            new table ( 158, 1.85327186949373479654360975305e280,  0 ),
            new table ( 159, 2.94670227249503832650433950735e282,  0 ),
            new table ( 160, 4.71472363599206132240694321176e284,  0 ),
            new table ( 161, 7.59070505394721872907517857094e286,  0 ),
            new table ( 162, 1.22969421873944943411017892849e289,  0 ),
            new table ( 163, 2.00440157654530257759959165344e291,  0 ),
            new table ( 164, 3.28721858553429622726333031164e293,  0 ),
            new table ( 165, 5.42391066613158877498449501421e295,  0 ),
            new table ( 166, 9.00369170577843736647426172359e297,  0 ),
            new table ( 167, 1.50361651486499904020120170784e300,  0 ),
            new table ( 168, 2.52607574497319838753801886917e302,  0 ),
            new table ( 169, 4.26906800900470527493925188890e304,  0 ),
            new table ( 170, 7.25741561530799896739672821113e306,  0 ),

            /*
            { 171, 1.24101807021766782342484052410e309,  0 },
            { 172, 2.13455108077438865629072570146e311,  0 },
            { 173, 3.69277336973969237538295546352e313,  0 },
            { 174, 6.42542566334706473316634250653e315,  0 },
            { 175, 1.12444949108573632830410993864e318,  0 },
            { 176, 1.97903110431089593781523349201e320,  0 },
            { 177, 3.50288505463028580993296328086e322,  0 },
            { 178, 6.23513539724190874168067463993e324,  0 },
            { 179, 1.11608923610630166476084076055e327,  0 },
            { 180, 2.00896062499134299656951336898e329,  0 },
            { 181, 3.63621873123433082379081919786e331,  0 },
            { 182, 6.61791809084648209929929094011e333,  0 },
            { 183, 1.21107901062490622417177024204e336,  0 },
            { 184, 2.22838537954982745247605724535e338,  0 },
            { 185, 4.12251295216718078708070590390e340,  0 },
            { 186, 7.66787409103095626397011298130e342,  0 },
            { 187, 1.43389245502278882136241112750e345,  0 },
            { 188, 2.69571781544284298416133291969e347,  0 },
            { 189, 5.09490667118697324006491921822e349,  0 },
            { 190, 9.68032267525524915612334651460e351,  0 },
            { 191, 1.84894163097375258881955918429e354,  0 },
            { 192, 3.54996793146960497053355363384e356,  0 },
            { 193, 6.85143810773633759312975851330e358,  0 },
            { 194, 1.32917899290084949306717315158e361,  0 },
            { 195, 2.59189903615665651148098764559e363,  0 },
            { 196, 5.08012211086704676250273578535e365,  0 },
            { 197, 1.00078405584080821221303894971e368,  0 },
            { 198, 1.98155243056480026018181712043e370,  0 },
            { 199, 3.94328933682395251776181606966e372,  0 },
            { 200, 7.88657867364790503552363213932e374,  0 }
            */
        };

        static table[] doub_fact_table = {
          new table ( 0,  1.000000000000000000000000000,    1L    ),
          new table ( 1,  1.000000000000000000000000000,    1L    ),
          new table ( 2,  2.000000000000000000000000000,    2L    ),
          new table ( 3,  3.000000000000000000000000000,    3L    ),
          new table ( 4,  8.000000000000000000000000000,    8L    ),
          new table ( 5,  15.00000000000000000000000000,    15L   ),
          new table ( 6,  48.00000000000000000000000000,    48L   ),
          new table ( 7,  105.0000000000000000000000000,    105L  ),
          new table ( 8,  384.0000000000000000000000000,    384L  ),
          new table ( 9,  945.0000000000000000000000000,    945L  ),
          new table ( 10, 3840.000000000000000000000000,    3840L   ),
          new table ( 11, 10395.00000000000000000000000,    10395L  ),
          new table ( 12, 46080.00000000000000000000000,       46080L       ),
          new table ( 13, 135135.0000000000000000000000,       135135L      ),
          new table ( 14, 645120.00000000000000000000000,      645120L      ),
          new table ( 15, 2.02702500000000000000000000000e6,   2027025L     ),
          new table ( 16, 1.03219200000000000000000000000e7,   10321920L    ),
          new table ( 17, 3.4459425000000000000000000000e7,    34459425L    ),
          new table ( 18, 1.85794560000000000000000000000e8,   185794560L   ),
          new table ( 19, 6.5472907500000000000000000000e8,            0 ),
          new table ( 20, 3.7158912000000000000000000000e9,            0 ),
          new table ( 21, 1.37493105750000000000000000000e10,          0 ),
          new table ( 22, 8.1749606400000000000000000000e10,           0 ),
          new table ( 23, 3.1623414322500000000000000000e11,           0 ),
          new table ( 24, 1.96199055360000000000000000000e12,          0 ),
          new table ( 25, 7.9058535806250000000000000000e12,           0 ),
          new table ( 26, 5.1011754393600000000000000000e13,           0 ),
          new table ( 27, 2.13458046676875000000000000000e14,          0 ),
          new table ( 28, 1.42832912302080000000000000000e15,          0 ),
          new table ( 29, 6.1902833536293750000000000000e15,           0 ),
          new table ( 30, 4.2849873690624000000000000000e16,           0 ),
          new table ( 31, 1.91898783962510625000000000000e17,          0 ),
          new table ( 32, 1.37119595809996800000000000000e18,          0 ),
          new table ( 33, 6.3326598707628506250000000000e18,           0 ),
          new table ( 34, 4.6620662575398912000000000000e19,           0 ),
          new table ( 35, 2.21643095476699771875000000000e20,          0 ),
          new table ( 36, 1.67834385271436083200000000000e21,          0 ),
          new table ( 37, 8.2007945326378915593750000000e21,           0 ),
          new table ( 38, 6.3777066403145711616000000000e22,           0 ),
          new table ( 39, 3.1983098677287777081562500000e23,           0 ),
          new table ( 40, 2.55108265612582846464000000000e24,          0 ),
          new table ( 41, 1.31130704576879886034406250000e25,          0 ),
          new table ( 42, 1.07145471557284795514880000000e26,          0 ),
          new table ( 43, 5.6386202968058350994794687500e26,           0 ),
          new table ( 44, 4.7144007485205310026547200000e27,           0 ),
          new table ( 45, 2.53737913356262579476576093750e28,          0 ),
          new table ( 46, 2.16862434431944426122117120000e29,          0 ),
          new table ( 47, 1.19256819277443412353990764062e30,          0 ),
          new table ( 48, 1.04093968527333324538616217600e31,          0 ),
          new table ( 49, 5.8435841445947272053455474391e31,           0 ),
          new table ( 50, 5.2046984263666662269308108800e32,           0 ),
          new table ( 51, 2.98022791374331087472622919392e33,          0 ),
          new table ( 52, 2.70644318171066643800402165760e34,          0 ),
          new table ( 53, 1.57952079428395476360490147278e35,          0 ),
          new table ( 54, 1.46147931812375987652217169510e36,          0 ),
          new table ( 55, 8.6873643685617511998269581003e36,           0 ),
          new table ( 56, 8.1842841814930553085241614926e37,           0 ),
          new table ( 57, 4.9517976900801981839013661172e38,           0 ),
          new table ( 58, 4.7468848252659720789440136657e39,           0 ),
          new table ( 59, 2.92156063714731692850180600912e40,       0 ),
          new table ( 60, 2.84813089515958324736640819942e41,       0 ),
          new table ( 61, 1.78215198865986332638610166557e42,       0 ),
          new table ( 62, 1.76584115499894161336717308364e43,       0 ),
          new table ( 63, 1.12275575285571389562324404931e44,       0 ),
          new table ( 64, 1.13013833919932263255499077353e45,       0 ),
          new table ( 65, 7.2979123935621403215510863205e45,        0 ),
          new table ( 66, 7.4589130387155293748629391053e46,        0 ),
          new table ( 67, 4.8896013036866340154392278347e47,        0 ),
          new table ( 68, 5.0720608663265599749067985916e48,        0 ),
          new table ( 69, 3.3738248995437774706530672060e49,        0 ),
          new table ( 70, 3.5504426064285919824347590141e50,        0 ),
          new table ( 71, 2.39541567867608200416367771623e51,       0 ),
          new table ( 72, 2.55631867662858622735302649017e52,       0 ),
          new table ( 73, 1.74865344543353986303948473285e53,       0 ),
          new table ( 74, 1.89167582070515380824123960272e54,       0 ),
          new table ( 75, 1.31149008407515489727961354964e55,       0 ),
          new table ( 76, 1.43767362373591689426334209807e56,       0 ),
          new table ( 77, 1.00984736473786927090530243322e57,       0 ),
          new table ( 78, 1.12138542651401517752540683649e58,       0 ),
          new table ( 79, 7.9777941814291672401518892225e58,        0 ),
          new table ( 80, 8.9710834121121214202032546920e59,        0 ),
          new table ( 81, 6.4620132869576254645230302702e60,        0 ),
          new table ( 82, 7.3562883979319395645666688474e61,        0 ),
          new table ( 83, 5.3634710281748291355541151243e62,        0 ),
          new table ( 84, 6.1792822542628292342360018318e63,        0 ),
          new table ( 85, 4.5589503739486047652209978556e64,        0 ),
          new table ( 86, 5.3141827386660331414429615754e65,        0 ),
          new table ( 87, 3.9662868253352861457422681344e66,        0 ),
          new table ( 88, 4.6764808100261091644698061863e67,        0 ),
          new table ( 89, 3.5299952745484046697106186396e68,        0 ),
          new table ( 90, 4.2088327290234982480228255677e69,        0 ),
          new table ( 91, 3.2122956998390482494366629620e70,        0 ),
          new table ( 92, 3.8721261107016183881809995223e71,        0 ),
          new table ( 93, 2.98743500085031487197609655470e72,       0 ),
          new table ( 94, 3.6397985440595212848901395509e73,        0 ),
          new table ( 95, 2.83806325080779912837729172696e74,       0 ),
          new table ( 96, 3.4942066022971404334945339689e75,        0 ),
          new table ( 97, 2.75292135328356515452597297515e76,       0 ),
          new table ( 98, 3.4243224702511976248246432895e77,        0 ),
          new table ( 99, 2.72539213975072950298071324540e78,       0 ),
          new table ( 100, 3.4243224702511976248246432895e79,       0 ),
          new table ( 101, 2.75264606114823679801052037785e80,      0 ),
          new table ( 102, 3.4928089196562215773211361553e81,       0 ),
          new table ( 103, 2.83522544298268390195083598919e82,      0 ),
          new table ( 104, 3.6325212764424704404139816015e83,       0 ),
          new table ( 105, 2.97698671513181809704837778865e84,      0 ),
          new table ( 106, 3.8504725530290186668388204976e85,       0 ),
          new table ( 107, 3.1853757851910453638417642339e86,       0 ),
          new table ( 108, 4.1585103572713401601859261374e87,       0 ),
          new table ( 109, 3.4720596058582394465875230149e88,       0 ),
          new table ( 110, 4.5743613929984741762045187512e89,       0 ),
          new table ( 111, 3.8539861625026457857121505465e90,       0 ),
          new table ( 112, 5.1232847601582910773490610013e91,       0 ),
          new table ( 113, 4.3550043636279897378547301176e92,       0 ),
          new table ( 114, 5.8405446265804518281779295415e93,       0 ),
          new table ( 115, 5.0082550181721881985329396352e94,       0 ),
          new table ( 116, 6.7750317668333241206863982681e95,       0 ),
          new table ( 117, 5.8596583712614601922835393732e96,       0 ),
          new table ( 118, 7.9945374848633224624099499564e97,       0 ),
          new table ( 119, 6.9729934618011376288174118541e98,       0 ),
          new table ( 120, 9.5934449818359869548919399477e99,       0 ),
          new table ( 121, 8.4373220887793765308690683435e100,      0 ),
          new table ( 122, 1.17040028778399040849681667362e102,       0 ),
          new table ( 123, 1.03779061691986331329689540625e103,       0 ),
          new table ( 124, 1.45129635685214810653605267528e104,       0 ),
          new table ( 125, 1.29723827114982914162111925781e105,       0 ),
          new table ( 126, 1.82863340963370661423542637086e106,       0 ),
          new table ( 127, 1.64749260436028300985882145742e107,       0 ),
          new table ( 128, 2.34065076433114446622134575470e108,       0 ),
          new table ( 129, 2.12526545962476508271787968008e109,       0 ),
          new table ( 130, 3.04284599363048780608774948111e110,       0 ),
          new table ( 131, 2.78409775210844225836042238090e111,       0 ),
          new table ( 132, 4.0165567115922439040358293151e112,        0 ),
          new table ( 133, 3.7028500103042282036193617666e113,        0 ),
          new table ( 134, 5.3821859935336068314080112822e114,        0 ),
          new table ( 135, 4.9988475139107080748861383849e115,        0 ),
          new table ( 136, 7.3197729512057052907148953438e116,        0 ),
          new table ( 137, 6.8484210940576700625940095873e117,        0 ),
          new table ( 138, 1.01012866726638733011865555744e119,       0 ),
          new table ( 139, 9.5193053207401613870056733264e119,        0 ),
          new table ( 140, 1.41418013417294226216611778042e121,       0 ),
          new table ( 141, 1.34222205022436275556779993902e122,       0 ),
          new table ( 142, 2.00813579052557801227588724819e123,       0 ),
          new table ( 143, 1.91937753182083874046195391280e124,       0 ),
          new table ( 144, 2.89171553835683233767727763739e125,       0 ),
          new table ( 145, 2.78309742114021617366983317355e126,       0 ),
          new table ( 146, 4.2219046860009752130088253506e127,        0 ),
          new table ( 147, 4.0911532090761177752946547651e128,        0 ),
          new table ( 148, 6.2484189352814433152530615189e129,        0 ),
          new table ( 149, 6.0958182815234154851890356000e130,        0 ),
          new table ( 150, 9.3726284029221649728795922783e131,        0 ),
          new table ( 151, 9.2046856051003573826354437561e132,        0 ),
          new table ( 152, 1.42463951724416907587769802630e134,       0 ),
          new table ( 153, 1.40831689758035467954322289468e135,       0 ),
          new table ( 154, 2.19394485655602037685165496051e136,       0 ),
          new table ( 155, 2.18289119124954975329199548675e137,       0 ),
          new table ( 156, 3.4225539762273917878885817384e138,        0 ),
          new table ( 157, 3.4271391702617931126684329142e139,        0 ),
          new table ( 158, 5.4076352824392790248639591467e140,        0 ),
          new table ( 159, 5.4491512807162510491428083336e141,        0 ),
          new table ( 160, 8.6522164519028464397823346347e142,        0 ),
          new table ( 161, 8.7731335619531641891199214170e143,        0 ),
          new table ( 162, 1.40165906520826112324473821082e145,       0 ),
          new table ( 163, 1.43002077059836576282654719098e146,       0 ),
          new table ( 164, 2.29872086694154824212137066574e147,       0 ),
          new table ( 165, 2.35953427148730350866380286512e148,       0 ),
          new table ( 166, 3.8158766391229700819214753051e149,        0 ),
          new table ( 167, 3.9404222333837968594685507847e150,        0 ),
          new table ( 168, 6.4106727537265897376280785126e151,        0 ),
          new table ( 169, 6.6593135744186166925018508262e152,        0 ),
          new table ( 170, 1.08981436813352025539677334714e154,       0 ),
          new table ( 171, 1.13874262122558345441781649128e155,       0 ),
          new table ( 172, 1.87448071318965483928245015709e156,       0 ),
          new table ( 173, 1.97002473472025937614282252992e157,       0 ),
          new table ( 174, 3.2615964409499994203514632733e158,        0 ),
          new table ( 175, 3.4475432857604539082499394274e159,        0 ),
          new table ( 176, 5.7404097360719989798185753611e160,        0 ),
          new table ( 177, 6.1021516157960034176023927864e161,        0 ),
          new table ( 178, 1.02179293302081581840770641427e163,       0 ),
          new table ( 179, 1.09228513922748461175082830877e164,       0 ),
          new table ( 180, 1.83922727943746847313387154568e165,       0 ),
          new table ( 181, 1.97703610200174714726899923887e166,       0 ),
          new table ( 182, 3.3473936485761926211036462131e167,        0 ),
          new table ( 183, 3.6179760666631972795022686071e168,        0 ),
          new table ( 184, 6.1592043133801944228307090322e169,        0 ),
          new table ( 185, 6.6932557233269149670791969232e170,        0 ),
          new table ( 186, 1.14561200228871616264651187999e172,       0 ),
          new table ( 187, 1.25163882026213309884380982464e173,       0 ),
          new table ( 188, 2.15375056430278638577544233437e174,       0 ),
          new table ( 189, 2.36559737029543155681480056857e175,       0 ),
          new table ( 190, 4.0921260721752941329733404353e176,        0 ),
          new table ( 191, 4.5182909772642742735162690860e177,        0 ),
          new table ( 192, 7.8568820585765647353088136358e178,        0 ),
          new table ( 193, 8.7203015861200493478863993359e179,        0 ),
          new table ( 194, 1.52423511936385355864990984535e181,       0 ),
          new table ( 195, 1.70045880929340962283784787050e182,       0 ),
          new table ( 196, 2.98750083395315297495382329688e183,       0 ),
          new table ( 197, 3.3499038543080169569905603049e184,        0 ),
          new table ( 198, 5.9152516512272428904085701278e185,        0 ),
          new table ( 199, 6.6663086700729537444112150067e186,        0 ),
          new table ( 200, 1.18305033024544857808171402556e188,       0 ),
          new table ( 201, 1.33992804268466370262665421635e189,       0 ),
          new table ( 202, 2.38976166709580612772506233164e190,       0 ),
          new table ( 203, 2.72005392664986731633210805920e191,       0 ),
          new table ( 204, 4.8751138008754445005591271565e192,        0 ),
          new table ( 205, 5.5761105496322279984808215214e193,        0 ),
          new table ( 206, 1.00427344298034156711518019425e195,       0 ),
          new table ( 207, 1.15425488377387119568553005492e196,       0 ),
          new table ( 208, 2.08888876139911045959957480403e197,       0 ),
          new table ( 209, 2.41239270708739079898275781478e198,       0 ),
          new table ( 210, 4.3866663989381319651591070885e199,        0 ),
          new table ( 211, 5.0901486119543945858536189892e200,        0 ),
          new table ( 212, 9.2997327657488397661373070276e201,        0 ),
          new table ( 213, 1.08420165434628604678682084470e203,       0 ),
          new table ( 214, 1.99014281187025170995338370390e204,       0 ),
          new table ( 215, 2.33103355684451500059166481610e205,       0 ),
          new table ( 216, 4.2987084736397436934993088004e206,        0 ),
          new table ( 217, 5.0583428183525975512839126509e207,        0 ),
          new table ( 218, 9.3711844725346412518284931849e208,        0 ),
          new table ( 219, 1.10777707721921886373117687056e210,       0 ),
          new table ( 220, 2.06166058395762107540226850068e211,       0 ),
          new table ( 221, 2.44818734065447368884590088393e212,       0 ),
          new table ( 222, 4.5768864963859187873930360715e213,        0 ),
          new table ( 223, 5.4594577696594763261263589712e214,        0 ),
          new table ( 224, 1.02522257519044580837604008002e216,       0 ),
          new table ( 225, 1.22837799817338217337843076851e217,       0 ),
          new table ( 226, 2.31700301993040752692985058084e218,       0 ),
          new table ( 227, 2.78841805585357753356903784452e219,       0 ),
          new table ( 228, 5.2827668854413291614000593243e220,        0 ),
          new table ( 229, 6.3854773479046925518730966640e221,        0 ),
          new table ( 230, 1.21503638365150570712201364459e223,       0 ),
          new table ( 231, 1.47504526736598397948268532937e224,       0 ),
          new table ( 232, 2.81888441007149324052307165546e225,       0 ),
          new table ( 233, 3.4368554729627426721946568174e226,        0 ),
          new table ( 234, 6.5961895195672941828239876738e227,        0 ),
          new table ( 235, 8.0766103614624452796574435210e228,        0 ),
          new table ( 236, 1.55670072661788142714646109101e230,       0 ),
          new table ( 237, 1.91415665566659953127881411447e231,       0 ),
          new table ( 238, 3.7049477293505577966085773966e232,        0 ),
          new table ( 239, 4.5748344070431728797563657336e233,        0 ),
          new table ( 240, 8.8918745504413387118605857518e234,        0 ),
          new table ( 241, 1.10253509209740466402128414180e236,       0 ),
          new table ( 242, 2.15183364120680396827026175195e237,       0 ),
          new table ( 243, 2.67916027379669333357172046456e238,       0 ),
          new table ( 244, 5.2504740845446016825794386748e239,        0 ),
          new table ( 245, 6.5639426708018986672507151382e240,        0 ),
          new table ( 246, 1.29161662479797201391454191399e242,       0 ),
          new table ( 247, 1.62129383968806897081092663913e243,       0 ),
          new table ( 248, 3.2032092294989705945080639467e244,        0 ),
          new table ( 249, 4.0370216608232917373192073314e245,        0 ),
          new table ( 250, 8.0080230737474264862701598667e246,        0 ),
          new table ( 251, 1.01329243686664622606712104019e248,       0 ),
          new table ( 252, 2.01802181458435147454008028642e249,       0 ),
          new table ( 253, 2.56362986527261495194981623168e250,       0 ),
          new table ( 254, 5.1257754090442527453318039275e251,        0 ),
          new table ( 255, 6.5372561564451681274720313908e252,        0 ),
          new table ( 256, 1.31219850471532870280494180544e254,       0 ),
          new table ( 257, 1.68007483220640820876031206743e255,       0 ),
          new table ( 258, 3.3854721421655480532367498580e256,        0 ),
          new table ( 259, 4.3513938154145972606892082546e257,        0 ),
          new table ( 260, 8.8022275696304249384155496309e258,        0 ),
          new table ( 261, 1.13571378582320988503988335446e260,       0 ),
          new table ( 262, 2.30618362324317133386487400329e261,       0 ),
          new table ( 263, 2.98692725671504199765489322224e262,       0 ),
          new table ( 264, 6.0883247653619723214032673687e263,        0 ),
          new table ( 265, 7.9153572302948612937854670389e264,        0 ),
          new table ( 266, 1.61949438758628463749326912007e266,       0 ),
          new table ( 267, 2.11340038048872796544071969939e267,       0 ),
          new table ( 268, 4.3402449587312428284819612418e268,        0 ),
          new table ( 269, 5.6850470235146782270355359914e269,        0 ),
          new table ( 270, 1.17186613885743556369012953528e271,       0 ),
          new table ( 271, 1.54064774337247779952663025366e272,       0 ),
          new table ( 272, 3.1874758976922247332371523360e273,        0 ),
          new table ( 273, 4.2059683394068643927077005925e274,        0 ),
          new table ( 274, 8.7336839596766957690697974006e275,        0 ),
          new table ( 275, 1.15664129333688770799461766294e277,       0 ),
          new table ( 276, 2.41049677287076803226326408256e278,       0 ),
          new table ( 277, 3.2038963825431789511450909263e279,        0 ),
          new table ( 278, 6.7011810285807351296918741495e280,        0 ),
          new table ( 279, 8.9388709072954692736948036845e281,        0 ),
          new table ( 280, 1.87633068800260583631372476186e283,       0 ),
          new table ( 281, 2.51182272495002686590823983534e284,       0 ),
          new table ( 282, 5.2912525401673484584047038284e285,        0 ),
          new table ( 283, 7.1084583116085760305203187340e286,        0 ),
          new table ( 284, 1.50271572140752696218693588728e288,       0 ),
          new table ( 285, 2.02591061880844416869829083919e289,       0 ),
          new table ( 286, 4.2977669632255271118546366376e290,        0 ),
          new table ( 287, 5.8143634759802347641640947085e291,        0 ),
          new table ( 288, 1.23775688540895180821413535163e293,       0 ),
          new table ( 289, 1.68035104455828784684342337075e294,       0 ),
          new table ( 290, 3.5894949676859602438209925197e295,        0 ),
          new table ( 291, 4.8898215396646176343143620089e296,        0 ),
          new table ( 292, 1.04813253056430039119572981576e298,       0 ),
          new table ( 293, 1.43271771112173296685410806860e299,       0 ),
          new table ( 294, 3.08150963985904315011544565835e300,       0 ),
          new table ( 295, 4.2265172478091122522196188024e301,        0 ),
          new table ( 296, 9.1212685339827677243417191487e302,        0 ),
          new table ( 297, 1.25527562259930633890922678431e304,       0 ),
          /*
          { 298, 2.71813802312686478185383230631e305,       0 },
          { 299, 3.7532741115719259533385880851e306,        0 },
          { 300, 8.1544140693805943455614969189e307,  }
          */
        };

        /* Chebyshev coefficients for Gamma*(3/4(t+1)+1/2), -1<t<1
        */
        static double[] gstar_a_data = {
          2.16786447866463034423060819465,
         -0.05533249018745584258035832802,
          0.01800392431460719960888319748,
         -0.00580919269468937714480019814,
          0.00186523689488400339978881560,
         -0.00059746524113955531852595159,
          0.00019125169907783353925426722,
         -0.00006124996546944685735909697,
          0.00001963889633130842586440945,
         -6.3067741254637180272515795142e-06,
          2.0288698405861392526872789863e-06,
         -6.5384896660838465981983750582e-07,
          2.1108698058908865476480734911e-07,
         -6.8260714912274941677892994580e-08,
          2.2108560875880560555583978510e-08,
         -7.1710331930255456643627187187e-09,
          2.3290892983985406754602564745e-09,
         -7.5740371598505586754890405359e-10,
          2.4658267222594334398525312084e-10,
         -8.0362243171659883803428749516e-11,
          2.6215616826341594653521346229e-11,
         -8.5596155025948750540420068109e-12,
          2.7970831499487963614315315444e-12,
         -9.1471771211886202805502562414e-13,
          2.9934720198063397094916415927e-13,
         -9.8026575909753445931073620469e-14,
          3.2116773667767153777571410671e-14,
         -1.0518035333878147029650507254e-14,
          3.4144405720185253938994854173e-15,
         -1.0115153943081187052322643819e-15
        };
        static ChebyshevSeries gstar_a_cs = new ChebyshevSeries(gstar_a_data, 29, -1, 1, 17);


        /* Chebyshev coefficients for
         * x^2(Gamma*(x) - 1 - 1/(12x)), x = 4(t+1)+2, -1 < t < 1
         */
        static double[] gstar_b_data = {
          0.0057502277273114339831606096782,
          0.0004496689534965685038254147807,
         -0.0001672763153188717308905047405,
          0.0000615137014913154794776670946,
         -0.0000223726551711525016380862195,
          8.0507405356647954540694800545e-06,
         -2.8671077107583395569766746448e-06,
          1.0106727053742747568362254106e-06,
         -3.5265558477595061262310873482e-07,
          1.2179216046419401193247254591e-07,
         -4.1619640180795366971160162267e-08,
          1.4066283500795206892487241294e-08,
         -4.6982570380537099016106141654e-09,
          1.5491248664620612686423108936e-09,
         -5.0340936319394885789686867772e-10,
          1.6084448673736032249959475006e-10,
         -5.0349733196835456497619787559e-11,
          1.5357154939762136997591808461e-11,
         -4.5233809655775649997667176224e-12,
          1.2664429179254447281068538964e-12,
         -3.2648287937449326771785041692e-13,
          7.1528272726086133795579071407e-14,
         -9.4831735252566034505739531258e-15,
         -2.3124001991413207293120906691e-15,
          2.8406613277170391482590129474e-15,
         -1.7245370321618816421281770927e-15,
          8.6507923128671112154695006592e-16,
         -3.9506563665427555895391869919e-16,
          1.6779342132074761078792361165e-16,
         -6.0483153034414765129837716260e-17
        };
        static ChebyshevSeries gstar_b_cs = new ChebyshevSeries(gstar_b_data, 29, -1, 1, 18);

        /* coefficients for gamma=7, kmax=8  Lanczos method */
        static double[] lanczos_7_c = {
          0.99999999999980993227684700473478,
          676.520368121885098567009190444019,
         -1259.13921672240287047156078755283,
          771.3234287776530788486528258894,
         -176.61502916214059906584551354,
          12.507343278686904814458936853,
         -0.13857109526572011689554707,
          9.984369578019570859563e-6,
          1.50563273514931155834e-7
        };

        /* Chebyshev expansion for System.Math.Log(gamma(x)/gamma(8))
         * 5 < x < 10
         * -1 < t < 1
         */
        static double[] gamma_5_10_data = {
         -1.5285594096661578881275075214,
          4.8259152300595906319768555035,
          0.2277712320977614992970601978,
         -0.0138867665685617873604917300,
          0.0012704876495201082588139723,
         -0.0001393841240254993658962470,
          0.0000169709242992322702260663,
         -2.2108528820210580075775889168e-06,
          3.0196602854202309805163918716e-07,
         -4.2705675000079118380587357358e-08,
          6.2026423818051402794663551945e-09,
         -9.1993973208880910416311405656e-10,
          1.3875551258028145778301211638e-10,
         -2.1218861491906788718519522978e-11,
          3.2821736040381439555133562600e-12,
         -5.1260001009953791220611135264e-13,
          8.0713532554874636696982146610e-14,
         -1.2798522376569209083811628061e-14,
          2.0417711600852502310258808643e-15,
         -3.2745239502992355776882614137e-16,
          5.2759418422036579482120897453e-17,
         -8.5354147151695233960425725513e-18,
          1.3858639703888078291599886143e-18,
         -2.2574398807738626571560124396e-19
        };
        static ChebyshevSeries gamma_5_10_cs = new ChebyshevSeries(gamma_5_10_data, 23, -1, 1, 11);
        #endregion

        #region static methods

        /* Lanczos method for real x > 0;
         * gamma=7, truncated at 1/(z+8) 
         * [J. SIAM Numer. Anal, Ser. B, 1 (1964) 86]
         */
        public static int LnGammaLanczos(double x, ref SpecialFunctionResult result)
        {
            int k;
            double Ag;
            double term1, term2;

            x -= 1.0; /* Lanczos writes z! instead of Gamma(z) */

            Ag = lanczos_7_c[0];
            for (k = 1; k <= 8; k++) { Ag += lanczos_7_c[k] / (x + k); }

            /* (x+0.5)*System.Math.Log(x+7.5) - (x+7.5) + LogRootTwoPi_ + System.Math.Log(Ag(x)) */
            term1 = (x + 0.5) * System.Math.Log((x + 7.5) / PZMath.M_E);
            term2 = LogRootTwoPi_ + System.Math.Log(Ag);
            result.Val = term1 + (term2 - 7.0);
            result.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * (System.Math.Abs(term1) + System.Math.Abs(term2) + 7.0);
            result.Err += PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val);

            return PZMath_errno.PZMath_SUCCESS;
        }

        /* x = eps near zero
         * gives double-precision for |eps| < 0.02
         */
        public static int LnGammaSgn0(double eps, ref SpecialFunctionResult lng, ref double sgn)
        {
            /* calculate series for g(eps) = Gamma(eps) eps - 1/(1+eps) - eps/2 */
            double c1 = -0.07721566490153286061;
            double c2 = -0.01094400467202744461;
            double c3 = 0.09252092391911371098;
            double c4 = -0.01827191316559981266;
            double c5 = 0.01800493109685479790;
            double c6 = -0.00685088537872380685;
            double c7 = 0.00399823955756846603;
            double c8 = -0.00189430621687107802;
            double c9 = 0.00097473237804513221;
            double c10 = -0.00048434392722255893;
            double g6 = c6 + eps * (c7 + eps * (c8 + eps * (c9 + eps * c10)));
            double g = eps * (c1 + eps * (c2 + eps * (c3 + eps * (c4 + eps * (c5 + eps * g6)))));

            /* calculate Gamma(eps) eps, a positive quantity */
            double gee = g + 1.0 / (1.0 + eps) + 0.5 * eps;

            lng.Val = System.Math.Log(gee / System.Math.Abs(eps));
            lng.Err = 4.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(lng.Val);
            sgn = System.Math.Sign(eps);

            return PZMath_errno.PZMath_SUCCESS;
        }


        /* x near a negative integer
         * Calculates sign as well as System.Math.Log(|gamma(x)|).
         * x = -N + eps
         * assumes N >= 1
         */
        public static int LnGammaSgnSing(int N, double eps, ref SpecialFunctionResult lng, ref double sgn)
        {
            if (eps == 0.0)
            {
                lng.Val = 0.0;
                lng.Err = 0.0;
                sgn = 0.0;
                PZMath_errno.ERROR("error", PZMath_errno.PZMath_EDOM);
            }
            else if (N == 1)
            {
                /* calculate series for
                 * g = eps gamma(-1+eps) + 1 + eps/2 (1+3eps)/(1-eps^2)
                 * double-precision for |eps| < 0.02
                 */
                double c0 = 0.07721566490153286061;
                double c1 = 0.08815966957356030521;
                double c2 = -0.00436125434555340577;
                double c3 = 0.01391065882004640689;
                double c4 = -0.00409427227680839100;
                double c5 = 0.00275661310191541584;
                double c6 = -0.00124162645565305019;
                double c7 = 0.00065267976121802783;
                double c8 = -0.00032205261682710437;
                double c9 = 0.00016229131039545456;
                double g5 = c5 + eps * (c6 + eps * (c7 + eps * (c8 + eps * c9)));
                double g = eps * (c0 + eps * (c1 + eps * (c2 + eps * (c3 + eps * (c4 + eps * g5)))));

                /* calculate eps gamma(-1+eps), a negative quantity */
                double gam_e = g - 1.0 - 0.5 * eps * (1.0 + 3.0 * eps) / (1.0 - eps * eps);

                lng.Val = System.Math.Log(System.Math.Abs(gam_e) / System.Math.Abs(eps));
                lng.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(lng.Val);
                sgn = (eps > 0.0 ? -1.0 : 1.0);
                return PZMath_errno.PZMath_SUCCESS;
            }
            else
            {
                double g;

                /* series for System.Math.Sin(Pi(N+1-eps))/(Pi eps) modulo the sign
                 * double-precision for |eps| < 0.02
                 */
                double cs1 = -1.6449340668482264365;
                double cs2 = 0.8117424252833536436;
                double cs3 = -0.1907518241220842137;
                double cs4 = 0.0261478478176548005;
                double cs5 = -0.0023460810354558236;
                double e2 = eps * eps;
                double sin_ser = 1.0 + e2 * (cs1 + e2 * (cs2 + e2 * (cs3 + e2 * (cs4 + e2 * cs5))));

                /* calculate series for ln(gamma(1+N-eps))
                 * double-precision for |eps| < 0.02
                 */
                double aeps = System.Math.Abs(eps);
                double c1, c2, c3, c4, c5, c6, c7;
                double lng_ser;
                SpecialFunctionResult c0 = new SpecialFunctionResult();
                SpecialFunctionResult psi_0 = new SpecialFunctionResult();
                SpecialFunctionResult psi_1 = new SpecialFunctionResult();
                SpecialFunctionResult psi_2 = new SpecialFunctionResult();
                SpecialFunctionResult psi_3 = new SpecialFunctionResult();
                SpecialFunctionResult psi_4 = new SpecialFunctionResult();
                SpecialFunctionResult psi_5 = new SpecialFunctionResult();
                SpecialFunctionResult psi_6 = new SpecialFunctionResult();
                psi_2.Val = 0.0;
                psi_3.Val = 0.0;
                psi_4.Val = 0.0;
                psi_5.Val = 0.0;
                psi_6.Val = 0.0;
                LnFactE(N, ref c0);
                Psi.PsiIntE(N + 1, ref psi_0);
                Psi.Psi1IntE(N + 1, ref psi_1);
                if (aeps > 0.00001) Psi.PsiNE(2, N + 1.0, ref psi_2);
                if (aeps > 0.0002) Psi.PsiNE(3, N + 1.0, ref psi_3);
                if (aeps > 0.001) Psi.PsiNE(4, N + 1.0, ref psi_4);
                if (aeps > 0.005) Psi.PsiNE(5, N + 1.0, ref psi_5);
                if (aeps > 0.01) Psi.PsiNE(6, N + 1.0, ref psi_6);
                c1 = psi_0.Val;
                c2 = psi_1.Val / 2.0;
                c3 = psi_2.Val / 6.0;
                c4 = psi_3.Val / 24.0;
                c5 = psi_4.Val / 120.0;
                c6 = psi_5.Val / 720.0;
                c7 = psi_6.Val / 5040.0;
                lng_ser = c0.Val - eps * (c1 - eps * (c2 - eps * (c3 - eps * (c4 - eps * (c5 - eps * (c6 - eps * c7))))));

                /* calculate
                 * g = ln(|eps gamma(-N+eps)|)
                 *   = -ln(gamma(1+N-eps)) + ln(|eps Pi/System.Math.Sin(Pi(N+1+eps))|)
                 */
                g = -lng_ser - System.Math.Log(sin_ser);

                lng.Val = g - System.Math.Log(System.Math.Abs(eps));
                lng.Err = c0.Err + 2.0 * PZMath_machine.PZMath_DBL_EPSILON * (System.Math.Abs(g) + System.Math.Abs(lng.Val));

                sgn = (PZMath.IsOdd(N) ? -1.0 : 1.0) * (eps > 0.0 ? 1.0 : -1.0);

                return PZMath_errno.PZMath_SUCCESS;
            }
            return PZMath_errno.PZMath_SUCCESS;
        }


        public static int LnGamma1Pade(double eps, ref SpecialFunctionResult result)
        {
            /* Use (2,2) Pade for Log[Gamma[1+eps]]/eps
             * plus a correction series.
             */
            double n1 = -1.0017419282349508699871138440;
            double n2 = 1.7364839209922879823280541733;
            double d1 = 1.2433006018858751556055436011;
            double d2 = 5.0456274100274010152489597514;
            double num = (eps + n1) * (eps + n2);
            double den = (eps + d1) * (eps + d2);
            double pade = 2.0816265188662692474880210318 * num / den;
            double c0 = 0.004785324257581753;
            double c1 = -0.01192457083645441;
            double c2 = 0.01931961413960498;
            double c3 = -0.02594027398725020;
            double c4 = 0.03141928755021455;
            double eps5 = eps * eps * eps * eps * eps;
            double corr = eps5 * (c0 + eps * (c1 + eps * (c2 + eps * (c3 + c4 * eps))));
            result.Val = eps * (pade + corr);
            result.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val);
            return PZMath_errno.PZMath_SUCCESS;
        }

        public static int LnGamma2Pade(double eps, ref SpecialFunctionResult result)
        {
            /* Use (2,2) Pade for Log[Gamma[2+eps]]/eps
             * plus a correction series.
             */
            double n1 = 1.000895834786669227164446568;
            double n2 = 4.209376735287755081642901277;
            double d1 = 2.618851904903217274682578255;
            double d2 = 10.85766559900983515322922936;
            double num = (eps + n1) * (eps + n2);
            double den = (eps + d1) * (eps + d2);
            double pade = 2.85337998765781918463568869 * num / den;
            double c0 = 0.0001139406357036744;
            double c1 = -0.0001365435269792533;
            double c2 = 0.0001067287169183665;
            double c3 = -0.0000693271800931282;
            double c4 = 0.0000407220927867950;
            double eps5 = eps * eps * eps * eps * eps;
            double corr = eps5 * (c0 + eps * (c1 + eps * (c2 + eps * (c3 + c4 * eps))));
            result.Val = eps * (pade + corr);
            result.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val);
            return PZMath_errno.PZMath_SUCCESS;
        }


        /* series for gammastar(x)
         * double-precision for x > 10.0
         */
        public static int GammaStarSer(double x, ref SpecialFunctionResult result)
        {
            /* Use the Stirling series for the correction to Log(Gamma(x)),
             * which is better behaved and easier to compute than the
             * regular Stirling series for Gamma(x). 
             */
            double y = 1.0 / (x * x);
            double c0 = 1.0 / 12.0;
            double c1 = -1.0 / 360.0;
            double c2 = 1.0 / 1260.0;
            double c3 = -1.0 / 1680.0;
            double c4 = 1.0 / 1188.0;
            double c5 = -691.0 / 360360.0;
            double c6 = 1.0 / 156.0;
            double c7 = -3617.0 / 122400.0;
            double ser = c0 + y * (c1 + y * (c2 + y * (c3 + y * (c4 + y * (c5 + y * (c6 + y * c7))))));
            result.Val = System.Math.Exp(ser / x);
            result.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * result.Val * System.Math.Max(1.0, ser / x);
            return PZMath_errno.PZMath_SUCCESS;
        }


        /* gamma(x) for x >= 1/2
         * assumes x >= 1/2
         */
        public static int GammaXGTHalf(double x, ref SpecialFunctionResult result)
        {
            if (x == 0.5)
            {
                result.Val = 1.77245385090551602729817;
                result.Err = PZMath_machine.PZMath_DBL_EPSILON * result.Val;
                return PZMath_errno.PZMath_SUCCESS;
            }
            else if (x <= (GSL_SF_FACT_NMAX + 1.0) && x == System.Math.Floor(x))
            {
                int n = (int)System.Math.Floor(x);
                result.Val = fact_table[n - 1]._f;
                result.Err = PZMath_machine.PZMath_DBL_EPSILON * result.Val;
                return PZMath_errno.PZMath_SUCCESS;
            }
            else if (System.Math.Abs(x - 1.0) < 0.01)
            {
                /* Use series for Gamma[1+eps] - 1/(1+eps).
                 */
                double eps = x - 1.0;
                double c1 = 0.4227843350984671394;
                double c2 = -0.01094400467202744461;
                double c3 = 0.09252092391911371098;
                double c4 = -0.018271913165599812664;
                double c5 = 0.018004931096854797895;
                double c6 = -0.006850885378723806846;
                double c7 = 0.003998239557568466030;
                result.Val = 1.0 / x + eps * (c1 + eps * (c2 + eps * (c3 + eps * (c4 + eps * (c5 + eps * (c6 + eps * c7))))));
                result.Err = PZMath_machine.PZMath_DBL_EPSILON;
                return PZMath_errno.PZMath_SUCCESS;
            }
            else if (System.Math.Abs(x - 2.0) < 0.01)
            {
                /* Use series for Gamma[1 + eps].
                 */
                double eps = x - 2.0;
                double c1 = 0.4227843350984671394;
                double c2 = 0.4118403304264396948;
                double c3 = 0.08157691924708626638;
                double c4 = 0.07424901075351389832;
                double c5 = -0.00026698206874501476832;
                double c6 = 0.011154045718130991049;
                double c7 = -0.002852645821155340816;
                double c8 = 0.0021039333406973880085;
                result.Val = 1.0 + eps * (c1 + eps * (c2 + eps * (c3 + eps * (c4 + eps * (c5 + eps * (c6 + eps * (c7 + eps * c8)))))));
                result.Err = PZMath_machine.PZMath_DBL_EPSILON;
                return PZMath_errno.PZMath_SUCCESS;
            }
            else if (x < 5.0)
            {
                /* Exponentiating the logarithm is fine, as
                 * long as the exponential is not so large
                 * that it greatly amplifies the error.
                 */
                SpecialFunctionResult lg = new SpecialFunctionResult();
                LnGammaLanczos(x, ref lg);
                result.Val = System.Math.Exp(lg.Val);
                result.Err = result.Val * (lg.Err + 2.0 * PZMath_machine.PZMath_DBL_EPSILON);
                return PZMath_errno.PZMath_SUCCESS;
            }
            else if (x < 10.0)
            {
                /* This is a sticky area. The logarithm
                 * is too large and the gammastar series
                 * is not good.
                 */
                double gamma_8 = 5040.0;
                double t = (2.0 * x - 15.0) / 5.0;
                SpecialFunctionResult c = new SpecialFunctionResult();
                ChebyshevSeries.ChebEvalE(gamma_5_10_cs, t, ref c);
                result.Val = System.Math.Exp(c.Val) * gamma_8;
                result.Err = result.Val * c.Err;
                result.Err += 2.0 * PZMath_machine.PZMath_DBL_EPSILON * result.Val;
                return PZMath_errno.PZMath_SUCCESS;
            }
            else if (x < GSL_SF_GAMMA_XMAX)
            {
                /* We do not want to exponentiate the logarithm
                 * if x is large because of the inevitable
                 * inflation of the error. So we carefully
                 * use System.Math.Pow() and System.Math.Exp() with exact quantities.
                 */
                double p = System.Math.Pow(x, 0.5 * x);
                double e = System.Math.Exp(-x);
                double q = (p * e) * p;
                double pre = PZMath.M_SQRT2 * PZMath.M_SQRTPI * q / System.Math.Sqrt(x);
                SpecialFunctionResult gstar = new SpecialFunctionResult();
                int stat_gs = GammaStarSer(x, ref gstar);
                result.Val = pre * gstar.Val;
                result.Err = (x + 2.5) * PZMath_machine.PZMath_DBL_EPSILON * result.Val;
                return stat_gs;
            }
            else
            {
                PZMath_errno.OverFlowError(ref result);
            }
            return PZMath_errno.PZMath_SUCCESS;
        }


        /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/


        public static int LnGammaE(double x, ref SpecialFunctionResult result)
        {

            if (System.Math.Abs(x - 1.0) < 0.01)
            {
                /* Note that we must amplify the errors
                 * from the Pade evaluations because of
                 * the way we must pass the argument, i.e.
                 * writing (1-x) is a loss of precision
                 * when x is near 1.
                 */
                int stat = LnGamma1Pade(x - 1.0, ref result);
                result.Err *= 1.0 / (PZMath_machine.PZMath_DBL_EPSILON + System.Math.Abs(x - 1.0));
                return stat;
            }
            else if (System.Math.Abs(x - 2.0) < 0.01)
            {
                int stat = LnGamma2Pade(x - 2.0, ref result);
                result.Err *= 1.0 / (PZMath_machine.PZMath_DBL_EPSILON + System.Math.Abs(x - 2.0));
                return stat;
            }
            else if (x >= 0.5)
            {
                return LnGammaLanczos(x, ref result);
            }
            else if (x == 0.0)
            {
                PZMath_errno.DomainError(ref result);
            }
            else if (System.Math.Abs(x) < 0.02)
            {
                double sgn = 0.0;
                return LnGammaSgn0(x, ref result, ref sgn);
            }
            else if (x > -0.5 / (PZMath_machine.PZMath_DBL_EPSILON * PZMath.M_PI))
            {
                /* Try to extract a fractional
                 * part from x.
                 */
                double z = 1.0 - x;
                double s = System.Math.Sin(PZMath.M_PI * z);
                double abss = System.Math.Abs(s);
                if (s == 0.0)
                {
                    PZMath_errno.DomainError(ref result);
                }
                else if (abss < PZMath.M_PI * 0.015)
                {
                    /* x is near a negative integer, -N */
                    if (x < int.MinValue + 2.0)
                    {
                        result.Val = 0.0;
                        result.Err = 0.0;
                        PZMath_errno.ERROR("error", PZMath_errno.PZMath_EROUND);
                    }
                    else
                    {
                        int N = -(int)(x - 0.5);
                        double eps = x + N;
                        double sgn = 0.0;
                        return LnGammaSgnSing(N, eps, ref result, ref sgn);
                    }
                }
                else
                {
                    SpecialFunctionResult lg_z = new SpecialFunctionResult();
                    LnGammaLanczos(z, ref lg_z);
                    result.Val = PZMath.M_LNPI - (System.Math.Log(abss) + lg_z.Val);
                    result.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val) + lg_z.Err;
                    return PZMath_errno.PZMath_SUCCESS;
                }
            }
            else
            {
                /* |x| was too large to extract any fractional part */
                result.Val = 0.0;
                result.Err = 0.0;
                PZMath_errno.ERROR("error", PZMath_errno.PZMath_EROUND);
            }
            return PZMath_errno.PZMath_SUCCESS;
        }


        public static int LnGammaSgnE(double x, ref SpecialFunctionResult result_lg, ref double sgn)
        {
            if (System.Math.Abs(x - 1.0) < 0.01)
            {
                int stat = LnGamma1Pade(x - 1.0, ref result_lg);
                result_lg.Err *= 1.0 / (PZMath_machine.PZMath_DBL_EPSILON + System.Math.Abs(x - 1.0));
                sgn = 1.0;
                return stat;
            }
            else if (System.Math.Abs(x - 2.0) < 0.01)
            {
                int stat = LnGamma2Pade(x - 2.0, ref result_lg);
                result_lg.Err *= 1.0 / (PZMath_machine.PZMath_DBL_EPSILON + System.Math.Abs(x - 2.0));
                sgn = 1.0;
                return stat;
            }
            else if (x >= 0.5)
            {
                sgn = 1.0;
                return LnGammaLanczos(x, ref result_lg);
            }
            else if (x == 0.0)
            {
                sgn = 0.0;

                PZMath_errno.DomainError(ref result_lg);
            }
            else if (System.Math.Abs(x) < 0.02)
            {
                return LnGammaSgn0(x, ref result_lg, ref sgn);
            }
            else if (x > -0.5 / (PZMath_machine.PZMath_DBL_EPSILON * PZMath.M_PI))
            {
                /* Try to extract a fractional
                  * part from x.
                  */
                double z = 1.0 - x;
                double s = System.Math.Sin(PZMath.M_PI * x);
                double abss = System.Math.Abs(s);
                if (s == 0.0)
                {
                    sgn = 0.0;
                    PZMath_errno.DomainError(ref result_lg);
                }
                else if (abss < PZMath.M_PI * 0.015)
                {
                    /* x is near a negative integer, -N */
                    if (x < int.MinValue + 2.0)
                    {
                        result_lg.Val = 0.0;
                        result_lg.Err = 0.0;
                        sgn = 0.0;
                        PZMath_errno.ERROR("error", PZMath_errno.PZMath_EROUND);
                    }
                    else
                    {
                        int N = -(int)(x - 0.5);
                        double eps = x + N;
                        return LnGammaSgnSing(N, eps, ref result_lg, ref sgn);
                    }
                }
                else
                {
                    SpecialFunctionResult lg_z = new SpecialFunctionResult();
                    LnGammaLanczos(z, ref lg_z);
                    sgn = (s > 0.0 ? 1.0 : -1.0);
                    result_lg.Val = PZMath.M_LNPI - (System.Math.Log(abss) + lg_z.Val);
                    result_lg.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result_lg.Val) + lg_z.Err;
                    return PZMath_errno.PZMath_SUCCESS;
                }
            }
            else
            {
                /* |x| was too large to extract any fractional part */
                result_lg.Val = 0.0;
                result_lg.Err = 0.0;
                sgn = 0.0;
                PZMath_errno.ERROR("error", PZMath_errno.PZMath_EROUND);
            }
            return PZMath_errno.PZMath_SUCCESS;
        }


        public static int GammaE(double x, ref SpecialFunctionResult result)
        {
            if (x < 0.5)
            {
                int rint_x = (int)System.Math.Floor(x + 0.5);
                double f_x = x - rint_x;
                double sgn_gamma = (PZMath.IsEven(rint_x) ? 1.0 : -1.0);
                double sin_term = sgn_gamma * System.Math.Sin(PZMath.M_PI * f_x) / PZMath.M_PI;

                if (sin_term == 0.0)
                {
                    PZMath_errno.DomainError(ref result);
                }
                else if (x > -169.0)
                {
                    SpecialFunctionResult g = new SpecialFunctionResult();
                    GammaXGTHalf(1.0 - x, ref g);
                    if (System.Math.Abs(sin_term) * g.Val * PZMath_machine.PZMath_DBL_MIN < 1.0)
                    {
                        result.Val = 1.0 / (sin_term * g.Val);
                        result.Err = System.Math.Abs(g.Err / g.Val) * System.Math.Abs(result.Val);
                        result.Err += 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val);
                        return PZMath_errno.PZMath_SUCCESS;
                    }
                    else
                    {
                        PZMath_errno.UnderFlowError(ref result);
                    }
                }
                else
                {
                    /* It is hard to control it here.
                     * We can only exponentiate the
                     * logarithm and eat the loss of
                     * precision.
                     */
                    SpecialFunctionResult lng = new SpecialFunctionResult();
                    double sgn = 0.0;
                    int stat_lng = LnGammaSgnE(x, ref lng, ref sgn);
                    int stat_e = Exp.ExpMultErrE(lng.Val, lng.Err, sgn, 0.0, ref result);
                    return PZMath_errno.ErrorSelect2(stat_e, stat_lng);
                }
            }
            else
            {
                return GammaXGTHalf(x, ref result);
            }
            return PZMath_errno.PZMath_SUCCESS;
        }


        public static int GammaStarE(double x, ref SpecialFunctionResult result)
        {

            if (x <= 0.0)
            {
                PZMath_errno.DomainError(ref result);
            }
            else if (x < 0.5)
            {
                SpecialFunctionResult lg = new SpecialFunctionResult();
                int stat_lg = LnGammaE(x, ref lg);
                double lx = System.Math.Log(x);
                double c = 0.5 * (PZMath.M_LN2 + PZMath.M_LNPI);
                double lnr_val = lg.Val - (x - 0.5) * lx + x - c;
                double lnr_err = lg.Err + 2.0 * PZMath_machine.PZMath_DBL_EPSILON * ((x + 0.5) * System.Math.Abs(lx) + c);
                int stat_e = Exp.ExpErrE(lnr_val, lnr_err, ref result);
                return PZMath_errno.ErrorSelect2(stat_lg, stat_e);
            }
            else if (x < 2.0)
            {
                double t = 4.0 / 3.0 * (x - 0.5) - 1.0;
                return ChebyshevSeries.ChebEvalE(gstar_a_cs, t, ref result);
            }
            else if (x < 10.0)
            {
                double t = 0.25 * (x - 2.0) - 1.0;
                SpecialFunctionResult c = new SpecialFunctionResult();
                ChebyshevSeries.ChebEvalE(gstar_b_cs, t, ref c);
                result.Val = c.Val / (x * x) + 1.0 + 1.0 / (12.0 * x);
                result.Err = c.Err / (x * x);
                result.Err += 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val);
                return PZMath_errno.PZMath_SUCCESS;
            }
            else if (x < 1.0 / PZMath_machine.PZMath_ROOT4_DBL_EPSILON)
            {
                return GammaStarSer(x, ref result);
            }
            else if (x < 1.0 / PZMath_machine.PZMath_DBL_EPSILON)
            {
                /* Use Stirling formula for Gamma(x).
                 */
                double xi = 1.0 / x;
                result.Val = 1.0 + xi / 12.0 * (1.0 + xi / 24.0 * (1.0 - xi * (139.0 / 180.0 + 571.0 / 8640.0 * xi)));
                result.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val);
                return PZMath_errno.PZMath_SUCCESS;
            }
            else
            {
                result.Val = 1.0;
                result.Err = 1.0 / x;
                return PZMath_errno.PZMath_SUCCESS;
            }
            return PZMath_errno.PZMath_SUCCESS;
        }


        public static int GammaInvE(double x, ref SpecialFunctionResult result)
        {
            if (x <= 0.0 && x == System.Math.Floor(x))
            { /* negative integer */
                result.Val = 0.0;
                result.Err = 0.0;
                return PZMath_errno.PZMath_SUCCESS;
            }
            else if (x < 0.5)
            {
                SpecialFunctionResult lng = new SpecialFunctionResult();
                double sgn = 0.0;
                int stat_lng = LnGammaSgnE(x, ref lng, ref sgn);
                if (stat_lng == PZMath_errno.PZMath_EDOM)
                {
                    result.Val = 0.0;
                    result.Err = 0.0;
                    return PZMath_errno.PZMath_SUCCESS;
                }
                else if (stat_lng != PZMath_errno.PZMath_SUCCESS)
                {
                    result.Val = 0.0;
                    result.Err = 0.0;
                    return stat_lng;
                }
                else
                {
                    return Exp.ExpMultErrE(-lng.Val, lng.Err, sgn, 0.0, ref result);
                }
            }
            else
            {
                SpecialFunctionResult g = new SpecialFunctionResult();
                int stat_g = GammaXGTHalf(x, ref g);
                if (stat_g == PZMath_errno.PZMath_EOVRFLW)
                {
                    PZMath_errno.UnderFlowError(ref result);
                }
                else
                {
                    result.Val = 1.0 / g.Val;
                    result.Err = System.Math.Abs(g.Err / g.Val) * System.Math.Abs(result.Val);
                    result.Err += 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val);
                    PZMath_errno.CheckUnderFlow(result);
                    return PZMath_errno.PZMath_SUCCESS;
                }
            }
            return PZMath_errno.PZMath_SUCCESS;
        } // GammaInvE()

        public static int TaylorCoeffE(int n, double x, ref SpecialFunctionResult result)
        {
            if (x < 0.0 || n < 0)
            {
                PZMath_errno.DomainError(ref result);
            }
            else if (n == 0)
            {
                result.Val = 1.0;
                result.Err = 0.0;
                return PZMath_errno.PZMath_SUCCESS;
            }
            else if (n == 1)
            {
                result.Val = x;
                result.Err = 0.0;
                return PZMath_errno.PZMath_SUCCESS;
            }
            else if (x == 0.0)
            {
                result.Val = 0.0;
                result.Err = 0.0;
                return PZMath_errno.PZMath_SUCCESS;
            }
            else
            {
                double log2pi = PZMath.M_LNPI + PZMath.M_LN2;
                double ln_test = n * (System.Math.Log(x) + 1.0) + 1.0 - (n + 0.5) * System.Math.Log(n + 1.0) + 0.5 * log2pi;

                if (ln_test < PZMath_machine.PZMath_LOG_DBL_MIN + 1.0)
                {
                    PZMath_errno.UnderFlowError(ref result);
                }
                else if (ln_test > PZMath_machine.PZMath_LOG_DBL_MAX - 1.0)
                {
                    PZMath_errno.OverFlowError(ref result);
                }
                else
                {
                    double product = 1.0;
                    int k;
                    for (k = 1; k <= n; k++)
                    {
                        product *= (x / k);
                    }
                    result.Val = product;
                    result.Err = n * PZMath_machine.PZMath_DBL_EPSILON * product;
                    PZMath_errno.CheckUnderFlow(result);
                    return PZMath_errno.PZMath_SUCCESS;
                }
            }
            return PZMath_errno.PZMath_SUCCESS;
        } // TaylorCoeffE()


        public static int FactE(int n, ref SpecialFunctionResult result)
        {

            if (n < 18)
            {
                result.Val = fact_table[n]._f;
                result.Err = 0.0;
                return PZMath_errno.PZMath_SUCCESS;
            }
            else if (n <= GSL_SF_FACT_NMAX)
            {
                result.Val = fact_table[n]._f;
                result.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val);
                return PZMath_errno.PZMath_SUCCESS;
            }
            else
            {
                PZMath_errno.OverFlowError(ref result);
            }
            return PZMath_errno.PZMath_SUCCESS;
        } // FactE()


        public static int DoubleFactE(int n, ref SpecialFunctionResult result)
        {
            if (n < 26)
            {
                result.Val = doub_fact_table[n]._f;
                result.Err = 0.0;
                return PZMath_errno.PZMath_SUCCESS;
            }
            else if (n <= GSL_SF_DOUBLEFACT_NMAX)
            {
                result.Val = doub_fact_table[n]._f;
                result.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val);
                return PZMath_errno.PZMath_SUCCESS;
            }
            else
            {
                PZMath_errno.OverFlowError(ref result);
            }
            return PZMath_errno.PZMath_SUCCESS;
        } // DoubleFactE()


        public static int LnFactE(int n, ref SpecialFunctionResult result)
        {
            /* CHECK_POINTER(result) */

            if (n <= GSL_SF_FACT_NMAX)
            {
                result.Val = System.Math.Log(fact_table[n]._f);
                result.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val);
                return PZMath_errno.PZMath_SUCCESS;
            }
            else
            {
                LnGammaE(n + 1.0, ref result);
                return PZMath_errno.PZMath_SUCCESS;
            }
        } // LnFactE()


        public static int LnDoubleFactE(int n, ref SpecialFunctionResult result)
        {

            if (n <= GSL_SF_DOUBLEFACT_NMAX)
            {
                result.Val = System.Math.Log(doub_fact_table[n]._f);
                result.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val);
                return PZMath_errno.PZMath_SUCCESS;
            }
            else if (PZMath.IsOdd(n))
            {
                SpecialFunctionResult lg = new SpecialFunctionResult();
                LnGammaE(0.5 * (n + 2.0), ref lg);
                result.Val = 0.5 * (n + 1.0) * PZMath.M_LN2 - 0.5 * PZMath.M_LNPI + lg.Val;
                result.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val) + lg.Err;
                return PZMath_errno.PZMath_SUCCESS;
            }
            else
            {
                SpecialFunctionResult lg = new SpecialFunctionResult();
                LnGammaE(0.5 * n + 1.0, ref lg);
                result.Val = 0.5 * n * PZMath.M_LN2 + lg.Val;
                result.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val) + lg.Err;
                return PZMath_errno.PZMath_SUCCESS;
            }
        } // LnDoubleFactE()


        public static int LnChooseE(int n, int m, ref SpecialFunctionResult result)
        {
            /* CHECK_POINTER(result) */

            if (m > n)
            {
                PZMath_errno.DomainError(ref result);
            }
            else if (m == n || m == 0)
            {
                result.Val = 0.0;
                result.Err = 0.0;
                return PZMath_errno.PZMath_SUCCESS;
            }
            else
            {
                SpecialFunctionResult nf = new SpecialFunctionResult();
                SpecialFunctionResult mf = new SpecialFunctionResult();
                SpecialFunctionResult nmmf = new SpecialFunctionResult();
                if (m * 2 > n)
                    m = n - m;
                LnFactE(n, ref nf);
                LnFactE(m, ref mf);
                LnFactE(n - m, ref nmmf);
                result.Val = nf.Val - mf.Val - nmmf.Val;
                result.Err = nf.Err + mf.Err + nmmf.Err;
                result.Err += 2.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val);
                return PZMath_errno.PZMath_SUCCESS;
            }
            return PZMath_errno.PZMath_SUCCESS;
        } // LnChooseE()


        public static int ChooseE(int n, int m, ref SpecialFunctionResult result)
        {
            if (m > n)
            {
                PZMath_errno.DomainError(ref result);
            }
            else if (m == n || m == 0)
            {
                result.Val = 1.0;
                result.Err = 0.0;
                return PZMath_errno.PZMath_SUCCESS;
            }
            else if (n <= GSL_SF_FACT_NMAX)
            {
                result.Val = (fact_table[n]._f / fact_table[m]._f) / fact_table[n - m]._f;
                result.Err = 6.0 * PZMath_machine.PZMath_DBL_EPSILON * System.Math.Abs(result.Val);
                return PZMath_errno.PZMath_SUCCESS;
            }
            else
            {
                if (m * 2 < n)
                    m = n - m;

                if (n - m < 64)  /* compute product for a manageable number of terms */
                {
                    double prod = 1.0;
                    int k;

                    for (k = n; k >= m + 1; k--)
                    {
                        double tk = (double)k / (double)(k - m);
                        if (tk > PZMath_machine.PZMath_DBL_MAX / prod)
                        {
                            PZMath_errno.OverFlowError(ref result);
                        }
                        prod *= tk;
                    }
                    result.Val = prod;
                    result.Err = 2.0 * PZMath_machine.PZMath_DBL_EPSILON * prod * System.Math.Abs(n - m);
                    return PZMath_errno.PZMath_SUCCESS;
                }
                else
                {
                    SpecialFunctionResult lc = new SpecialFunctionResult();
                    int stat_lc = LnChooseE(n, m, ref lc);
                    int stat_e = Exp.ExpErrE(lc.Val, lc.Err, ref result);
                    return PZMath_errno.ErrorSelect2(stat_lc, stat_e);
                }
            }
            return PZMath_errno.PZMath_SUCCESS;
        } // ChooseE()

        /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

        public static double LnGamma(double x)
        {
            SpecialFunctionResult result = new SpecialFunctionResult();
            LnGammaE(x, ref result);
            return result.Val;
        } // LnGamma()
        #endregion
    }
}
